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Abstract 

We consider an inverse problem arising in thermo-/photo- acoustic tomography that amounts 
to reconstructing a function / from its circular or spherical means with the centers lying on a 
given measurement surface. (Equivalently, these means can be expressed through the solution 
u{t,x) of the wave equation with the initial pressure equal to /.) An explicit solution of this 
inverse problem is obtained in 3D for the surface that is the boundary of an open octet, and 
in 2D for the case when the centers of integration circles lie on two rays starting at the origin 
and intersecting at the angle equal to tt/N, N = 2,3,4,.... Our formulas reconstruct the Radon 
projections of a function closely related to /, from the values of u(t, x) on the measurement 
surface. Then, function / can be found by inverting the Radon transform. 

Keywords: photoacoustic tomography, thermoacoustic tomography, optoacoustic tomography, 
wave equation, spherical means, explicit inversion formulas 

1 Introduction 

Thermoacoustic tomography (TAT) [17,41] and photo- or opto- acoustic tomography (PAT or 
OAT) [3,16,34] are based on the thermoacoustic effect: when a material is heated it expands. 
To perform measurements, a biological object is illuminated with a short electromagnetic pulse 
that heats the tissue and, through the resulting thermoacoustic expansion, generates an outgoing 
acoustic wave. Acoustic pressure is then measured by detectors placed on a surface (completely or 
partially) surrounding the object. The source reconstruction inverse problem of TAT/PAT/OAT 
consists in finding the initial pressure within the object from the measurements. This pressure is 
closely related to the blood content in the tissues, making these imaging modalities particularly 
effective for cancer detection and imaging of blood vessels. 

Under the assumption that the speed of sound in the tissues is constant, for certain simple 
acquisition surfaces a solution of this inverse problem can be represented by explicit closed-form 
formulas. Such formulas are similar to the well-known filtration/backprojection formula in X-ray 
tomography; they allow one to compute the initial pressure by evaluating an explicit integro- 
differential operator at each node of a reconstruction grid. The existence and form of explicit 
inversion formulas are closely related to the shape of the data acquisition surface. For the simplest 
case of a planar surface, explicit formulas have been known for several decades [1, 7,9]. The authors 
of [42] found a filtration/backprojection formula that is valid for a plane, a 3D sphere and an 
infinite 3D cylinder. In [10,11, 21, 33] several different inversion formulas were derived for spherical 
acquisition surfaces in spaces of various dimensions. In [23] explicit reconstruction formulas were 
obtained for detectors placed on a surface of a cube (or square, in 2D) or on surfaces of certain other 
polygons and polyhedra. Several authors [14, 28, 35, 38] recently found explicit inversion formulas 
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for the data measured from surfaces of an ellipse (in 2D) or an ellipsoid (in 3D) surrounding 
the object. These formulas can be further extended by continuity to an elliptic paraboloid or 
parabolic cylinder [29, 30]. Certain more complicated polynomial surfaces (including a paraboloid) 
were considered in [35], although not all of these surfaces are attractive from a practical point of 
view, as they would require surrounding the object by several layers of detectors. In addition to 
the explicit formulas, there exist several reconstruction algorithms (for closed acquisition surfaces) 
based on various series expansions [13,22,24,31,32]. These techniques sometimes lead to very 
fast implementations (e.g. [22, 24]); however, their efficient numerical implementation may require 
certain non-trivial computational skills. 

In this paper we derive inversion formulas for acquisition schemes that use either linear arrays 
of point- or line- detectors, or two dimensional arrays of sensors partially surrounding the object. 
Several experimental setups motivate our research. For example, in [5], linear detectors were 
arranged in a corner-like assembly rotated along the object. Such a scheme requires inversion of 
series of 2D circular Radon transforms with centers located on two perpendicular rays. A similar 3D 
problem sometime arises when optically scanned planar glass plate is used as an acoustic detector, 
as done by researchers from University College London (see, for example, [6, 8]). The use of a single 
glass plate provides only a limited angle view of the object, leading to disappearance of some edges 
in the image. To eliminate such artifacts, one needs to repeat scanning with a detector repositioned 
to view the object at a different angle, or to add additional glass plates and/or acoustical mirrors. 
The methods of the present paper are applicable in the case of consecutive measurements. The 
acquisition scheme with multiple glass detectors may need to be modeled as a reverberating cavity 
due to multiple reflections of acoustic waves from the glass surfaces. The methods of the present 
paper are not applicable in such situation; instead, the reader is referred to [6, 8,15, 20]. 

Recently, linear arrays of point-like piezoelectric detectors started to gain popularity due to 
their fast rate of acquisition (see [18, 25, 26]). Similarly to planar arrays, one such linear assembly 
yields only limited-angle measurements, and one has to either use mirrors (as is done in [25]), or 
to arrange several arrays around the object (reflections from linear arrays can be neglected due to 
their smaller cross-section). 

In this paper we consider a 3D acquisition scheme where point detectors are located on the 
boundary of an infinite octant, and a similar 2D problem with detectors covering the boundary of 
an inhnite angular sector with an opening angle vr/A^, N = 2,3,4,.... In both cases it is assumed 
that reflection of acoustic waves from the detectors is either absent or negligible. Unlike many of 
inversion formulas mentioned above, the formulas we derive are not of a filtration/backprojection 
type. Instead, we obtain explicit expressions allowing one to recover from the measurements the 
classical Radon projections of the unknown initial pressure. The latter function then is easily 
reconstructed by application of the well-known inverse Radon transform. 

2 Preliminaries 

2.1 Formulation of the problem 

We assume throughout the paper that the attenuation of acoustic waves is negligible, and that the 
speed of sound in the tissues is constant (and, thus, it can be set to 1 without loss of generality). 
This simplified model is acceptable in such important applications as breast imaging, and most of 
the formulas and algorithms mentioned in the Introduction are based on these assumptions. Then, 
the acoustic pressure p{t, x) solves the following initial problem for the wave equation in the whole 
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Figure 1: Acquisition geometry in 3D and 2D 


space d = 3, 

( ptt = Ap, X G t £ [ 0 ,oo), 

< p{0,x) = f{x), ( 1 ) 

[ Pti0,x) = 0. 

were the initial pressure of the acoustic wave f{x) is the function we seek. It is assumed that f{x) 
is finitely supported within the region of interest D, and that the point-wise measurements of p{t, x) 
are done outside D. A very similar two-dimensional problem arises when linear integrating detectors 
are used for measurements (e.g., see [13]). For example, if a set of linear detectors parallel to vector 
63 = ( 0 , 0 , 1 ) is utilized, they measure integrals /(t,xi,X 2 ) of the pressure along lines parallel to 
63 and the problem consists in finding the initial distribution of these integrals /(x) = /(O, xi,X 2 ). 
This problem can also be formulated in the form (1), with d = 2, with p{t, x) replaced by I{t, xi, X 2 ). 

Our goal is to reconstruct /(x) from measurements of p{t, x) made at all points of the measuring 
surface (or, in 2D, measuring curve). In 3D, as the measuring surface we will use the boundary dQ 
of the first Cartesian octet Q in 

Q = {x £ M^jxi > 0, X 2 > 0, X 3 > 0}. ( 2 ) 

For the 2D problem the data u{t, x) will be assumed known on the boundary of infinite angular 
segment Q contained between the ray M"'' = {(c, 0)| c > 0} and ray = {(ccos/3, csin/3)| c > 0}, 
where angle j3 is to one of the values n/N, N = 2,3,4,... . In this case the measuring curve dQ is 
the union M'*' U Figure 1 shows schematically the measuring sets in 2D and 3D. In both cases 
it is assumed that the support D of / is properly contained in Q. 

Thus, our ultimate goal is to find explicit reconstruction formulas expressing /(x) through 
measured values p{t,y), y £ dQ, t £ [0,oo). 

2.2 Outline of the approach 

Solution of the initial value problem (1) can be written as a convolution of / with the time derivative 
of the retarded (or causal) Green’s function G+(t,x) of the wave equation in free space [40]: 

P{t,y) =-^ J f{x)G^{t,y - x)dx. (3) 

Q 
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In 3D the above-mentioned Green’s function takes form 


where 5{t) is the Dirac’s delta function. In 2D 


5{t — |a:|) 
47r|x| 



0 


t > |xl, 
otherwise 


The measurements p{t, y) are just the values of the convolution (3) at points y lying on dQ. The 
time anti-derivative of the data P{t,y) can be easily recovered from p{t,y): 


P{t, y) = ^(0’ 2/) = 0> 2/ € dQ. 


Function P{t,y) is related to f{x) by the convolution 


P{t,y) = J f{x)G^{t,y - x)dx, yGdQ. 
n 

For notational convenience we extend p{t, y) to negative values of t by zero: 

p(t,y) = 0, t < 0. 


(4) 


(With this convention equation (3) is valid for all t G M.) 

Below we will also use the advanced free-space Green’s function G~ {t, x) that allows one to 
solve wave equation backwards in time. It can be defined simply by reversing the sign in front of t: 

G~{t,x) = G~^{—t,x). 


Gonsider now an integrable function ip{t,y) defined on WxdQ, finitely supported in the first 
variable and decaying in y fast enough so that the integral below converges absolutely. One can 
define a retarded single layer potential uj {t, x) with density y) by the following equation [12, 39]; 


t-\y-x\ 


= J j p{T,y)G^{t-T,y-x)dTdy 


dQ -oo 

ip{t - s, y)G+(s, y - x)dsdy. 

dQ \y-x\ 

Similarly, an advanced single layer potential is given by the formula 

-\y-A 


Uu,it,x) = 


p{t-s,y)G {s,y-x)d.sdy 


dQ -oo 

OO 


dQ \y-x\ 


J ifitP s,y)G^{s,y - x)dsdy. 
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Due to vanishing of y — x) when s < \y — x\, one can equivalently write 


u^{t,x) = 


(fit =F s, 2 /)G+(s, y - x)dsdy. 


dQ K 


Our approach to the inverse problem in hand is based on the observation that the integral 
< f,u~{0, •) > defined as 

>= j fix)u~{0,x)dx 
n 

can be easily reconstructed from the data P{t,y): 


/ f{x)u,^{0,x)dx = 

[ fix) 

/ / ^{s,y)G^{s,y - x)dsdy 

J 

n 

J 

o 

dQ R 


J J P{s,y) 


dQ R 


f{x)G^{s, y — x)dx 


dx 


dsdy 


ip{s,y)P{s,y)dsdy. 


dQ R 


This procedure can be repeated for a set of single layer potentials u~(t,x) defined by different 
densities (p. If these potentials evaluated at t = 0 form some sort of a basis, the above integrals 
provide projections of / on the basis elements, making it possible to reconstruct /. 

Finding families of basis functions represented by single layer potentials is not necessarily easy, 
especially in unbounded domains. In general, certain solutions of the wave equation can be rep¬ 
resented by a combination of a single- and double- layer potentials. For example, if u{t,x) is a 
finitely supported in time solution of the wave equation in a bounded region Q* with a boundary 
dQ*, then u can be represented in the following form [12, 39] 


u{t, x) = 


A. 

dn 


i{t =F s, y)G+{s, y - x)dsdy 


dQ* R 


J J u{t^ s,y)-^G^{s,y - x)dsdy. 


dQ* R 


where second term is a double layer potential with density u{T,y). Sign =F in the above equation 
indicates that u(t, x) can be expressed through its boundary values (including the values of the 
normal derivative) either in the past, at the times (t — s), or in the future (at the time values t + s). 

However, for our purposes we need representations expressed only by a single potential sup¬ 
ported on a boundary of a domain, and our domains are not bounded. In what follows we show 
that such representations can be obtained for certain combinations of delta-like plane waves. This 
will allow us to recover from the data the integrals of products of f{x) with these functions; such 
integrals are related to the values of the Radon transform of f{x). Then, f(x) can be reconstructed 
by inverting the Radon transform. 
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3 Single layer representations 

In this section we find single layer representations for certain combinations of plane waves. 

3.1 2D case 

Consider now a Coxeter cross, i.e. a set of N straight lines passing through the origin, with the 
angle (3 = tt/N between any two adjacent lines; without loss of generality we assume that one 
of these lines coincides with xi-axis. The boundary of the 2D angular segment Q arising in the 
2D formulation of our inverse problem (see Section 2.1) is a subset of this set of points. (For a 
discussion of injectivity of the spherical mean Radon transform with centers lying on a Coxeter 
cross we refer the reader to [2].) 

3.1.1 An “odd” operator 

In this section we define an operator that maps an arbitrary continuous function of x G into a 
function that is odd with respect to reflections about all the lines constituting the Coxeter cross. 

Define a reflector 91 mapping x = (xi,X 2 )^ onto 91x = (xi,—X 2 )^, and a linear operator T^p 
rotating any vector x G by the angle (p counterclockwise. Obviously 

^2^ = = I (5) 

where I is the identity transformation. Now, let us define an operator O transforming a continuous 
function h{x) into a continuous function ho{x) according to the following formula 

N-l 

ho{x) = Oh{x) = Y, [K'^ 2 px) - hmlpx)] . ( 6 ) 

A:=0 

It follows from (5) and (6) that ho{x) is invariant with respect to rotations by the angle 2/3 

ho{^ 20 x) = ho{T~Y) = ho{x), 

and, in general 

ho (j^ 2 / 3 x'^ = ho{x), k eZ. ( 7 ) 

Further, we notice that due to (5) 

N-l N-l N-l 

Y hmlf,x) = Y hm-^x) = Y h{T%^x), ( 8 ) 

k=i) k=0 k=0 

where the second equality is verified as follows. Each x can be written as \x\T 0 e 1 , where ei = (1,0)^ 
and 9 is some angle. Then 

91T-^"x = 91T-j^ (|x|T,ei) = |x|91T(_2fc^+,)ei = |x|T( 2 fc^_,)ei = T|^91 x, 

which proofs the second equality in (8). Therefore, operator O (defined by equation (6)) admits an 
alternative representation 

N-l 

ho{x) = Oh{x) = Y [h{'^%x) - h{T%^x)\ . (9) 

A:=0 
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Since = I, it is immediately clear from (9) that 

hoi^x) = -ho{x), 

i.e., ho{x) is an odd function in X 2 - Moreover, due to the rotation invariance (7), ho{x) is odd 
with respect to any line that is a rotation of the xi axis by the angle 2/3k, k £ and 

hoi^T'^/sx) = -ho{x), k£Z. (10) 

One can also show that ho is odd with respect to any line that is a rotation of xi-axis by 
the angle {2k + l)/3, k £ "L. Let’s consider two points cT^T^ei and cT-qT^ci symmetric with 
respect to the line tT^ei, t £ (— 00 , 00 ). Due to periodicity (equation (5)), ho admits yet another 
representation 

7V-1 

ho{x) = Oh{x) = \h{Tlpx) - /i(T^+^inx) . 
k=0 

Then 

7V-1 

hoicT^T^ei) = Y [h{cTlf^xT^T0ei)-h{cY+^mT^Tpei) 

k=0 

7V-1 

“ ^ [H(^^2l3{k+l)+a^l) — h{(^^2/3{k+l)-a^l)] = “^o(cT_aT,gei). 
k=0 

Therefore, ho is odd with respect to a rotation of xi-axis by the angle (3, and due to the 2/3- 
periodicity (equation (7)), it is odd with respect to a rotation of xi-axis by the angle {2k + l)/3. 
Thus, we have proved the following 

Proposition 1 Operator O defined by equation (6) maps a eontinuous function into a continuous 
function that is odd with respect to reflections about each of the lines constituting the Coxeter eross. 

Corollary 2 If ho{x) = Oh{x), then ho{x) vanishes on eaeh of the lines eonstituting the Coxeter 
eross, and, in partieular, on the boundary dQ of segment Q. 

Proposition 3 If support of h{x) is eontained within Q, then restriction of ho to Q coincides with 
h, i.e. 

h{x) = ho{x), Vx £ Q. 

Proof. If support of h{x) is contained within Q, then each term in (6) that is not h{x) itself, is 
supported within one of rotated and/or reflected versions of segment Q not intersecting Q itself. ■ 

Proposition 4 Suppose g{x) and h{x) are eontinuous functions, one of them is eompaetly sup¬ 
ported, with go{x) = Og{x) and ho{x) = Oh{x). Then 

/ go{x)h{x)dx = / ho{x)g{x)dx. 

R2 R2 
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Proof. By using the definition of operator O and interchanging the order of summation and 
integration one obtains 


Af-l 


N-\ 


go{x)h{x)dx = ^ /” g{'^2gx)h{x)dx - ^ / g{^T20x)h{x)dx 

fc= 0 jg 2 ^=°R 2 

N-1 „ N-1 „ 

= X] / g{x)h{T~^x)dx - ^ / g{x)h{d\T~^x)dx 

1 _r\ 7,_n 


J g{x) I [h(T 


-k 

2/3 


x) - himrgx) 


dx = 


k=0 

N-1 


g{x) f ^ \dx= ho{x)g{x)dx, 


.k=0 


where (5) was used on the last line. ■ 

Finally, the following statement follows directly from (7) and (10). 

Proposition 5 If g{x) and h{x) are continuous functions and g{x) is compactly supported then 


/ g{x)ho dx= g{x)ho{x)dx, 

R2 R2 

/ g{x)ho dx = — g{x)ho (x) dx, /c G Z. 

R2 R2 


3.2 3D case 

In the 3D version of our problem, region Q is the first Cartesian octet Q defined by equation (2). 

3.2.1 An “odd” operator 

In this section we define an operator that maps an arbitrary continuous function of 3D variable 
into function that is odd with respect to reflections about any of the three coordinate planes. 

Let us define operators 9Ij, j = 1,2,3, reflecting any x = (xi,/C 2 ,a/ 3 )^ with respect to the 
coordinate plane with the normal ej: 

iRlX = (-Xi,X2,X3)^, 

9123/ = (xi,-X 2 ,X 3 )^, (11) 

9I3X = (xi,X2, -X3)^. 

It follows from the above dehnition that these reflectors commute. 

We can now define an operator O mapping continuous function h{x) into function ho{x) odd 
with respect to reflections about each of the three coordinate planes: 


ho{x) = Oh{x) = h{x) — h{d\ix) — h{d\ 2 x) — /i(9l3x) 

+ /i(9l291ix) + /i(9l391ix) + h{m3d\2x) - h{d\3^R2^ix). 


(12) 



Since operators d\j commute and are self-invertible, verifying that operator O have the above- 
mentioned property is straightforward. 

Proposition 6 Operator O defined by equation (12) maps a eontinuous function into a continuous 
function that is odd with respect to reflections about each of the three eoordinate planes. 

Similarly, one proves 

Proposition 7 If h{x) is bounded and ho{x) = Oh{x) then 

ho i^ix) = -ho{x), ho = ho{x), ho = -ho{x), i,j,k G {1,2,3}. 

Corollary 8 If ho{x) = Oh{x), then ho{x) vanishes on each of the coordinate planes, in particu¬ 
lar, on the boundary dQ of the first octet Q. 

Proposition 9 If support of h{x) is eontained within Q, then restriction of ho to Q coincides with 
h, i.e. 

h{x) = ho{x), Vx G Q. 

Proposition 10 Suppose g{x) and h{x) are continuous functions, one of them is compactly sup¬ 
ported, with go{x) = Og{x) and ho{x) = Oh{x). Then 

/ go{x)h{x)dx = / ho{x)g{x)dx. 

K3 R3 

The proof of the last proposition is similar to that of Proposition 4. 

As a direct consequence of Proposition 7, one obtains 


Proposition 11 If g{x) and h{x) are continuous functions, and g{x) is compactly supported, then 

[g{x)ho {^)dx = - [g{x)ho{x)dx, 


g{x)ho (iHjiHjx) = jg{x)ho{x)dx, 


g{x)ho = - g{x)ho{x)dx, i,j,k G {1,2,3}. 


3.3 Waves represented by single layer potentials 
3.3.1 Odd combinations of waves 

Consider a non-negative, infinitely smooth, even function g (t) defined on R, compactly supported 
with all its derivatives on (—1,1), and such that f^ipdt = 1. Define scaled function gs{t) by the 
formula 

„^(t) = h (h) . 
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Figure 2: Constructing U(t,x) from u^^oj(t,x) in the cases N = 2 and N = 3. The dashed line 
shows the lines of the Coxeter cross; U{t, x) is odd with respect to these lines 


Family of functions %(t) is delta-approximating, in the sense of distributions: 

limrye(t) = (5(t). (13) 

£^0 

Now, let us consider a plane wave Ue(t,x) dehned by the equation 

Ue,uj(t,x) = r]s{t - X ■ Uj) , 

where cj is a unit vector defining the direction of propagation of such a wave. It is easy to check 
that Ue{t,x) is a solution of the wave equation in the whole space 

^2 

X t G R, 

and that, in the limit e ^ 0, approximates a delta-wave 

lim Ueu}{t,x) = 6{t — X ■ uj). 

£^0 ’ 

Let us now use the operator O introduced in the previous sections and define the following 
combination of plane waves: 

Ue,uj{t,x) = OUe^Ut^x). (14) 

To simplify the notation, within this section we will suppress the subscripts and denote this function 
by U{t,x). Function U{t,x) consists of the plane wave U£^u}{t,x) together with its reflections and 
(in 2D case) rotations, selected in such a way that it is odd with respect to a set of lines (or planes) 
containing the boundary of dQ. Thus, 

U{t,x) = 0, Mx^dQ, Vt G R. 


From now on, let us assume that vector uj is pointing strictly inside Q. In the 2D case the 
operator O is defined by equation (6). 

Apply it to U£^i^{t,x): 


U{t,x) = Ou£^aj{t,x) = 

N-l 

k=0 

N-l 

k=0 


r]£ (^t - 


rjelt-x- 


(15) 
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It follows that U{t,x) has the following form 


U{t, x) 

= rje {t — X ■ uj) + U*(t, x), 

(16) 


2N-1 


U*{t,x) 

= ^ (t-X- , 

(17) 


j=i 


where aj equal 1 or —1, and unit vectors j = 1, 2N — 1, are rotated and/or reflected versions 
of u). By construction of operators iH and T^ 2/3 entering equation (15), each of vectors lies in 
its own sector of the Coxeter cross. Figure 2 shows two examples of constructing U{t,x) from 
x), for the cases N = 2 and N = 3. 

In the 3D case, function U{t,x) has a similar form, given by equations (16), (17) with = 4; 
each of the vectors in this case is pointing inside its own octet of the Cartesian grid. 


3.3.2 Representing odd waves by single layer potentials 

We would like to represent U{t,x) in Q by an advanced single layer potential supported on dQ. 

Let us first consider a bounded subset Q* of Q defined as follows. Let i?(0, R) be the open ball 
of radius R centered at the origin. R will be used as a large parameter. Define Q* = Q (1 B{0, R). 
Boundary dQ* of Q* consists of two parts: 


dQ* = dQl U dQ*2, 
dQl = dQnB{0,R) 
dQ* = QndB{0,R). 


Since is compactly supported on (—e,e), function U{t,x) identically vanishes in B{0,R) for all 
t such that t > to = R + s. Therefore, for all x in Q* 


U{t, x) 


dQ* \y-x\ 


G^{s,y- x)^U{t + s,y)dsdy - U{t + s,y)-^G^{s,y 


OO 

= J J G+{s,y - x)-^U{t + s,y)dsdy 

dQl \y-x\ 


+ 


dQl \y-x\ 


d d 

G^{s,y- x)—U{t + s,y)dsdy - U{t + s,y)—G^{s,y 



dsdy 



dsdy, 


(18) 


where second equality is due to vanishing of U{t,y) on dQ and, thus, on dQ^. Let us denote by I 
the double integral on the last line of the above equation. Then, due to (16) 


1 = 


dQl \y-x\ 

OO 


d d 

G+(s, y - x)—r]e {t + s - y ■ u) dsdy - rje {t + s - y ■ uj) —G+(s, y - x) 
dn dn 


dsdy 


+ 


dQl \y-x\ 


G+(s, y - x)^U*{t + s, y)dsdy -U*{t + s, y)^G+(s, y - x) 


dsdy. 
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Let us show that the last double integral in the above equation vanishes if R is chosen sufficiently 
large. Let us first consider only values of x lying within region A = Q r\ B{0,ro), and values of t 
such that \t\ < tq. We notice that for fixed x and y, integration in time variable is done only over 
values t + s, with s > \x — y\. Consider values of an arbitrary term in U* (say, with subscript j as 
defined by equation (17)) at the point y and at time t + s. It has form {t + s — y ■ . Function 

ye (t + s — y ■ and its normal derivative vanish outside of the interval 

-e + y- <t + s<£ + y- 

On the other hand, for x £ A, the smallest value of \y — x\ is R — xq, so that integration is done 
over the interval 

s > R — Xq 


or 

t + s > R — 2 xq. 

It is sufficient to show that there is R such that 

£ + y ■ < R — 2ro (19) 

for all y G dQ^- Any such y can be represented as y = ^R, where C is a unit vector lying strictly 
within sector Q. Since each j = 1, ...,2N — 1 is fixed and lies strictly outside Q, there is a 
uniform (in y and j) upper bound a on the dot product C, ■ 

<a<l. ( 20 ) 

Therefore, for the left hand side of (19) we obtain 

£ + y ■ < e + Ra. 


It is easy to check that if R is chosen so that 


R > 


£ + 2X0 
(1-a)’ 


( 21 ) 


then inequality (19) is satisfied for each y G dQ 2 , x £ A, and t such that \t\ < xq. Therefore, for 
such X and t equation (18) simplifies to 


OO 

U{t,x)= J J G^{s,y 

dQt \y-x\ 


x)—U{t + s,y)dsdy + 


OO 



dQl \y-x\ 


G+{s,y 


^ d 


(t + s 


- Ve {t + s 


yuj) -^G+{s,y 



dsdy. 


y-uj) 


( 22 ) 


Let us now show that the second double integral in (22) equals % (t + s — y • cu) for \x\ < xq and 
|t| < Xq. Consider again ball B{0, R) with the boundary 5(0, R). Within this ball xj^ (t + s — y ■ lo) 
can be represented by the integral 


Vs = 


OO 

/ 


d 

G^{s,y - x)—ye {t + s - y ■ uj) - yit + s-y 


<-^)-^G+{s,y-x) 


dsdy. (23) 


S{0,R) \y-x\ 
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Considerations similar to the above, show that if R is chosen so that 

£ + ^ 

(1-/3)’ 

with some /? < 1, then the inner integral in (23) vanishes for all values y = RC, such that 

c • W < /3. 


(24) 


If one selects /3 close enough to 1 to guarantee that all y G satisfy ( ■ uj < j3, and if R satisfies 
(24), then 


OO 


dQ* \y-x\ 


G+{s,y 


\ ^ I 


y • w) - ?7e (t -h s 


y.u,)-^G*{s,y 



dsdy. 


(25) 


Now, by comparing (22) and (25) one obtains 


U{t,x) - r]e{t - X ■ u)) = J j G^{s,y - x)-^U{t + s,y)dsdy, |x| < ro, |t| < ro, (26) 

dQl \y-x\ 

valid if R satisfies both (21) and (24). Finally, we notice that for t G [—ro, —e] support of the wave 
%{t — X ■ uj) does not intersect Q, and therefore this term in (26) can be dropped: 


OO 

U{t,x)= J J G^{s,y-x)-^U{t + s,y)dsdy, 


9QI \y-x\ 


x\ < ro, t G [-ro, -e], x gQ. 


Since for a fixed ui this representation is valid for any sufficiently large R, one can take a limit i? —)• oo 
and replace the integration over dQl by integration over dQ. Moreover, since ro is arbitrary, this 
representation is valid within any bounded region 17 properly contained in Q, for t < —e: 


OO 

U{t,x) = J J G~^{s,y — x)—U{t + s,y)dsdy, t <—e, xGQ. 


dQ \y-x\ 


(27) 


Thus, we have proven 


Proposition 12 Let U{t,x) be a combination of plane waves defined by equations (16) and (17), 
with infinitely smooth % supported on the interval [—e,e] (where e > 0 is arbitrary) and vanishing 
with its derivatives at the endpoints of the interval. Then U{t,x) can be represented within any 
bounded region Q properly contained in Q by an advanced single layer potential (27); the represen¬ 
tation is valid for all t < —e. 

Let us make one more observation on the single layer representation (27). Spatial integration 
in this formula can be restricted to a bounded subset of dQ, per the following 


Proposition 13 Let function U{t,x) be defined as in the previous proposition, with fixed function 
rjs, and fixed values e > 0 and tv G Q. For arbitrary value ro > e > 0 consider a bounded subset Q* = 
Q (1 B(0,ro) of Q. Then, there exists value R = i?(a;,e,ro) such that the following representation 
of U {t, x) holds within Q*: 

OO 

U{t,x)= J J G^{s,y - x)-^U{t +s,y)dsdy, -ro < t <-e, x G Q*. (28) 

dQnB{0,R)) \y-x\ 
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Proof. The proof of this statement is quite similar to the proof of the previous proposition. Here 
is a brief sketch of it. All the points y on the integration surface dQ can be represented in the form 
y = Cl?/|- Since vector uj and its reflections/rotations lie correspondingly strictly within Q and 
within reflections/rotations of Q, there is an upper bound a{uj), such that an inequality similar to 
(20) is satisfied, with R = \y\ and a = a{uj). Therefore, for values of y such that 

li'i P'’* 

the inner integral in (27) vanishes, since support of + s,y) is disjoint from the integration 

interval. Hence, equation (28) holds. ■ 

4 Reconstruction of Radon projections 

The goal of this section is to recover the Radon projections of / from measurements p{t,y). More 
precisely, we will be able to reconstruct projections of Of] this, however, is enough to reconstruct 
/ by inverting the Radon transform. 

The latter transform, TZ, of an arbitrary continuous compactly supported function h{x) can be 
defined as follows [27]; 

{Tlh){t,uj)= / h{x)S{t — X ■ uj)dx, d = 2,3, (30) 

where unit vector uj is lying on a unit circle (in 2D) or a sphere (in 3D). For a fixed values of uj 
and t the Radon transform equals to an integral of a function over a hyperplane normal to uo and 
lying at a distance \t\ from the origin. The range of variable t is chosen so that all hyperplanes 
intersecting the support of h are accounted for. Inversion formulas allowing one to reconstruct h 
from known values of IZh are well known, along with a number of efficient reconstruction algorithms 
(e.g., [19,27]). 

Following the idea presented in Section 2.2, we utilize the single layer representation (27) and 
find the integral of f{x) with the combination of plane waves Ue,ui{t,x) defined by equation (14); 


f{x)Ue,Ut,x)dx = J f{x) 



d 

G+(s, y - x)—Ue,to{t + s, y)dsdy 


dQ \y-x\ 


dx 


A 

dn 


dQ 0 


Ue,coit + S,y) 


d 


f{x)G~^{s, y — x)dx 


dsdy 


P{s, y)-^Ue,uj{t + s, y)dsdy. 


dQ 0 


Let us modify formulas (17) and (16) and write 

27V-1 

^jVe (t- X- UJ^P'^ , 

1=0 


(31) 
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with (To = 1 an(d = oj. Now the normal (derivative of can be computed explicitly: 

2N-1 


d 


dn{y) 


Ue,uj{t + s,y) = n{y) ■ VUe,uj{t + s,y) = - ^ aj (n{y) ■ v'e (t + s - y ■ 


j=0 


ye{t + s-y- . 


2N-1 


j=0 

Substitute the above expression in (31) and integrate by parts in s; 

r2iv-i 


n 


fix)Ue,uj{t,x)dx = J j p{s,y) 

dQ 0 


^ Gj [n{y) ■ ys(t + s-y- 


j=0 


dsdy. 


By taking the limit e —)■ 0 and recalling (4) and (13) we obtain 

2Af-l 


lim / f{x)Ue,uj{t,x)dx= / 
£->■0 / / 
n dQ 


^ Gj {n{y) ■ p{y ■ - t, y) 


j=0 


dy. 


The left hand side of the above equation can be modified using Proposition 4: 

lim / f{x)U£uj{t,x)dx = lim / f (x)Oui;{t,x)dx = lim / uJt,x)Of{x)dx 
£^0 J ’ £^>0 J £-)-0 J 

nun 

= lim / u^{t, x)Of{x)dx = lim / Of(x)r]s (t — x ■ oj) dx 
£^0 J £-)-0 J 


n 


= / Of{x)6{t — x-oj)dx = {TZ[fo])it,oj), t<0, oj^Q, 


where fo{x) = Of{x), which proves 

Proposition 14 Suppose f{x) is a continuous function compactly supported in C Q, vector oj 
lies strictly in Q, and t < 0. Then, Radon projections (Jl[fo]) of foix) can be reconstructed 

from the boundary values of solution p{-,y) of the problem (1) by the formula 


{n[fo]){t,uj) = 


dQ 


2N-1 


^ Gj {n{y) ■ ) p{y ■ oj^^^ - t, y) 

3=0 


,(i) 


dy, 


(32) 


where coefficients gj and vectors oj^^'l have the same values as in (16) and (17). 


Corollary 15 Suppose function f{x), coefficients gj, and vectors ij^d) and uj are the same as in 
the previous proposition. Suppose, in addition, that region Tt d Q is contained within a ball B(0,ro) 
of radius tq and centered at 0. Then there is a value R (given by equation (29) with s = 0) such 
that the integration in (32) can be restricted to a finite subset of the boundary dQ H B{0, R{cj)): 


dy, ro <t < 0. (33) 


r 

2N-1 

i'^[fo]){t,uj) = J 

^ Gj (n{y) ■ p{y ■ - t, y) 

dQnB(0,R{uj)) 

_ 3=0 
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In order to find projections for directions of u not lying within Q, we observe that the Radon 
transform of an odd function Of{x) has many redundancies. First, it is well known [27] and is 
easily seen from the definition (30) that Radon transform TZh of a general function h is redundant 
in that 

(7^/i) {-t, -u) = {nh) {t, u). (34) 

Let us show that the Radon transform of an odd function Of{x) has additional redundancies. 
In 2D, using Proposition 5, one obtains 

(^[/o]) = lim J f{x)U^ -^k^Jt,x)dx = lim J f{x)Oue,u^ dx 

= lim [ f{x)Oue,uj (t, x) dx = (7^ [fo]) (t, u) , (35) 

e^Oj 

n 


and, similarly 

(7^ [fo]) {t, = - (7^ [fo]) {t, iv) . (36) 

In 3D, with the help of Proposition 11, we observe the following symmetries 

r (7^[/o])(^,lHiu;) = -(7^[/o])(^,u;) 

(7^[/o])(^,fn^^H,•a;) = (7^[/o])(^,w) i,j,kG {1,2,3}. (37) 

[ (7^ [fo]) (t, = - (7^ [fo]) (t, cu) 

Therefore, after projections (7Z [fo]) (t,^) of the function fo{x) are reconstructed using Proposition 
14 from the data p{t, y) for all cu G Q and all t < 0, one can recover projections for other values of 
u using formulas (35) and (36) (in 2D) or formulas (37) (in 3D). Values of {Tl[fo]) (L*^) for t > 0 
are then recovered using equation (34). For directions of cu parallel to the lines of Coxeter cross (in 
2D) or to coordinate planes (in 3D), Radon transform {'R-[fo]) (L^) vanishes. This follows from 
the fact that fo is odd with respect to reflections about the lines of Coxeter cross (in 2D) or about 
the coordinate planes (in 3D). Therefore, integration over hyperplanes orthogonal to these lines or 
planes yields zero. Finally, since both sides of the equation (32) are well-defined and continuous at 
t = 0, this formula extends to t = 0 by continuity. 

Thus, we have proven the following 

Theorem 16 Suppose f{x) is a continuous function compactly supported within a proper subset D 
of Q. Then, the Radon projections {TZ [fo]) {t,a!), t ^ R, u ^ of the function fo{x) = Of{x) 

can be explicitly reconstructed from the boundary values of solution p{-,y) of the problem (1), using 
formulas (32)-(36) (in 2D), or formulas (32), (34) and (37) (in 3D). If uj is parallel to the lines 
of Coxeter cross (in 2D) or to coordinate planes (in 3D), then {IZ[fo]) {t,uj) = 0, t G R. 

Finally, once projections {TZ[fo]) (L^^) are found for all values of t and all values of u lying 
on the unit circle (in 2D) or on the unit sphere (in 3D), function foix) can be reconstructed 
by application of one of the many known inversion techniques for the Radon transform (see, for 
example [19,27]). As mentioned above, function f(x) we seek to recover is just a restriction of 
fo{x) to Q. 

In addition. Corollary 15 implies that if data p{t, y) are known only on a finite subset dQ n 
B{0,R{uj)) of the boundary of Q, and / is supported within a ball of radius ro centered at the 


16 



origin, projections {Tl[fo]) can be reconstructed for those directions u that point inside Q 
and satisfy the inequality R > 2ro/(l — a(cj)), or 


a{u) < 1 


2ro 

R ■ 


Parameter a{uj) in the above inequality was defined in Section 3.3.2 as the upper bound on the dot 
products of vectors with all unit vectors lying in dQ. Due to the symmetries in the set of 
j = 0, ..2N — 1, and compactness of the set dQ n 


aiuj) = max (cv ■ C). 
C=9Qn§‘*-i 


Therefore, for fixed R and ro, and for t > 0, projections (TZ [fo]) can be reconstructed in the 
interval of values of uo satisfying the inequality 


max (a; • C) < 1 
C=9Qn§‘i-i 


2ro 
R ' 


In 3D, this inequality describes a ’’triangle-like” connected region T{R,rQ) that is a subset of the 
set of all unit vectors with positive coordinates. In 2D, if vector oj is represented as (cosy,sin 7 ) 
with 7 G [0,/3], (recall that f3 = tt/N, N = 2,3,4,...), the above inequality defines the interval 
I{R,ro) = [ 70 ,/3 — 7 o], where 70 = arccos (1 — 2ro/i?). By using the symmetries expressed by 
equations (34)-(37), one then finds values of {TZ [fo]) that corresponds to rotations and/or 

reflections of I{R,ro) or T{R,ro), and to nonnegative values of t. 


Theorem 17 Suppose f{x) is a continuous function compactly supported within a proper subset 
D 0 / (5 n B(0,ro), and boundary values of solution p{-,y) of the problem (1) are known for all 
y G dQ* = dQ n B{0, R),where ro and R are arbitrary (but with R > 2ro). Then, in 3D the 
Radon projections {TZ[fo]) {t,co) of the function fo{x) = Of{x) can be explicitly reconstructed 
from p{-,y), y G dQ*, for all values of oj = (0:1,002,013) such that max \ojj\ < 1 — and all t 

with |t| < ro, using equations (33), (34) and (37). In 2D, the Radon projections {IZ[fo\)(t,oj(y)) 

of the function fo(x) = Of{x) can be explicitly reconstructed from p(-,y), y G dQ*, for all 7 G 
2N-1 

U [70 + — 7o + all t with |t| < rg, using formulas (33) and (34)-(36). (In all 

k=^ 

dimensions ilZ [fo]) it,oj) = 0 for t > ro.) 


5 Numerical simulations 

In this section our formulas for reconstructing Radon projections are tested in two numerical 
simulations, one in 2D and the other one in 3D. The goal of these numerical experiments is to 
illustrate the exactness of the formulas, rather than to model realistic thermo- or photo- acoustic 
reconstruction. To this end, we did not simulate measurement noise, and did not study the effects of 
truncation of the acquisition surfaces. Moreover, when modeling the forward problem, an effort was 
made to sample solution of the wave equation on a grid covering the whole infinite boundary dQ. 
Clearly, this cannot be done numerically on a uniform grid. We, thus, used grids that are uniform 
on parts of the boundary relatively close to the origin (e.g. |y| is less than several hundreds), and 
that are getting progressively coarser for larger values of |y|. 

As stated by Corollary 15, for a particular direction of 00 integration is done only over a finite 
subset of dQ. The remote parts of dQ are only needed for computing projections with 00 near 
parallel to one of the parts of the boundary. Since remote parts of dQ were sampled very coarsely, 
values of reconstructed projections corresponding to such directions of oo were less accurate. 
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(e) Reconstruction from p{t,y) with \y\ < 300 (f) Reconstruction from noisy data, \y\ < 300 


Figure 3: Example of reconstructing projections from p{t, y) in 2D, N = 3 

5.1 A corner-like domain in 2D 

Consider a problem of reconstructing Radon projections from solution of the wave equation known 
on dQ consisting of two rays with angle /3 = vr/S between them. This geometry corresponds to 
a subset of Coxeter cross with N equal to 3. As a phantom f{x), we use a linear combination 
f{x) = hi{x) — h 2 {x) of functions hj, j = 1, 2, in the form 

hj = h (^ 7 ^) > ( 38 ) 

where h{s) is a smooth function with a maximum equal to 1, at s = 0, and vanishing with several of 
its derivatives at s = 1. This phantom is plotted in figure 3(a). The odd part fo{x) of the function 
f{x) is defined by formula (6); it is shown in figure 3(b). 

Values of p{t,y) corresponding to the chosen phantom f{x) were computed on dQ by first 
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(d) Projections on ( — 1,0) x (0,7r/2), ip = ipo (e) Recon. error on (—1,0) x (0,7r/2), ip = po 


Figure 4: Reconstructing projections {TZfo) from values ofp{t,y) known on dQ in 3D 


computing circular averages M{y, r) defined by the equation 


M{y,r) = J 
§1 


_ f f{y + ra) 




da, 


and by exploiting a well known connection between p{t,y) and M{y,r) : 


P{t,y) 


^ J f{x)G^{t,y - x)dx = -^ J f{y + z)G^{t,z)dz 

O 12 


--f f 

2TTdt J J 

t §1 


f{y + ra) 

\/t‘^ - kp 


da rdr 


Id r M{y,r) 

t ^ 


The Abel transform of M{y,r) and its time derivative were computed numerically. 

For further comparison, the Radon projections {Tlifo]) i't,‘x{9)) of fo{x) were computed first 
using the exact /o(x).The result is shown in figure 3(c). One can clearly see in this figure symmetries 
described by equations (34)-(36). Next, we reconstructed the values of projections fromp(t, y) using 
formulas (32)-(36). On a gray scale (or color) plot the result is indistinguishable from the one in 
figure 3(c). The plot of the difference between the reconstructed and known projections is shown 
in figure 3(d). One can notice that the error for values of 6 close to A: G Z, is larger than for 
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the rest of these values. An explanation was presented right before the start of this section. In an 
Loo-norm, the relative error over the whole range of values t and 9 is about 3%. If values of 9 close 
to ^ are excluded from the comparison, this error drops to about 1%. At least a part of this error 
is due to the error in computing y) by a three-step process involving finding circular averages, 
computing the Abel transform, and subsequent numerical differentiation. 

In figure 3(e) we present an image reconstructed from data p{t,y) known only for values \y\ < 
300. According to Theorem 17 applied with values ro = 1, i? = 300, our formulas yield exact 

5 

reconstruction only for values of cj( 7 ) with 7 G (J [6.62°-|-/c60°, (A:-|-1)60° — 6.62°]. If one examines 

k=0 

the image in part (e) only within this set of values, the error is similar to that shown in figure 3(d). 
For other values of 7 the error is larger, and can be noticed in a high-resolution version of the 
figure. 

Figure 3(f) demonstrates reconstruction from noisy (and truncated in y) data. A 100% noise (in 
L 2 norm) was added to values p{t, y) at each point y with \y\ < 300, on the time intervals of length 
2 starting at the time of arrival of the acoustic waves to point y. The effect of noise is noticeable 
in the figure; the relative Loo norm of the error is about 25%. This error, however, is quite small 
taking into account the high level of noise in the data — as one would expect in view of the highly 
stable nature of our reconstruction formulas. 

5.2 Boundary of an octet in 3D 

In the second numerical experiment we reconstructed Radon projections of a function from solution 
of the wave equation known on a boundary dQ of the first octet in 3D. As a phantom f(x) defined 
within the intersection of the unit ball and the first octet, we used a linear combination of four 
functions in the form (38). The centers xj of these functions lied either in the plane z = 0.3 or plane 
z = 0.6. The cross-sections of the phantom by these two planes are shown in figure 4, parts (a) 
and (b). Projections (77 [/o]) of function foix) with w parametrized in spherical coordinates 
(i.e., u) = {simp cos 9, sijnp sin 9, cos ip)) are shown in figure 4(c) for a fixed value of 99 = (^0 ~ 53°. 
Symmetries expressed by equations (34) and (37) can be clearly seen in this figure. A magnified 
image of projections restricted to the region {t,9) G [—1,0] x [0, 7 r/ 2 ] with (^ = (^0 is shown in 
figure 4(d). 

Modeling the forward problem (i.e., computing values of p{t, y) on dQ) is actually simpler in 3D 
than in 2D. For radial functions in the form (38) there exists an explicit exact formula [13] yielding 
values of p{t, y). We precomputed these values on polar grids defined on each of three infinite faces 
constituting dQ. The discretization step in the radial direction was uniform for relatively small 
values of the radius but became increasingly coarser for larger values. 

Numerical reconstruction of projections from p{t,y) performed using formulas (32), (34) and 
(37) in the region {t,9) G [—1,0] x [0, 7 r/ 2 ] yields data that are very close to the exact ones for 
most values of cp. The results obtained for ip = ipQ are indistinguishable from the exact values on 
a gray-scale or color plot. The error for ip = ipQ is plotted in figure 4(e); the relative error in Loo 
norm is about 0.1%. Similar error is observed for wide range of values of ip except those close to 0 
or 7 r/ 2 . We explain a better accuracy obtained in 3D simulation by a better accuracy in modeling 
p{t,y): in the 3D case these values are given by a theoretically exact formula. 

6 Concluding remarks and acknowledgement 

We presented explicit formulas that reconstruct Radon projections of fo from the values of the 
wave equation known on a boundary of an angular domain with vr/iV, N = 2,3,4,.. (in 2D) or on 
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the boundary of octet Q (in 3D). Here are some additional thoughts: 

• Since the solution of the wave equation in the case of constant speed of sound can be expressed 
through the spherical/circular means (and vice versa), finding Radon projections from these 
means is straightforward. 

• Explicit relations between the Radon transform and the spherical/circular means were found 
in [4] in the case when the measurement surface is a hyperplane. The present results may be 
considered a generalization of that work, although the case of a hyperplane in 2D or 3D is 
not presented in our paper. 

• General approach of this paper, as well as that of [4], is somewhat similar to the approach of 
Popov and Sushko [36, 37] who proposed a formula for finding Radon projections of a function 
from optoacoustic measurements. Their formula, however, yields a parametrix of the problem 
rather than the exact solution. 

• Theorem 16 can be easily generalized to the case when the acquisition surface (in 3D) is the 
boundary of an infinite region bounded by three planes, with two of these planes intersecting 
at the angle tt/N, N = 3,4, 5... and the third plane perpendicular to the first two. 

Acknowledgement: The author’s research is partially supported by the NSF, award NSF/DMS- 
1211521. 
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